#include "Smoothing.h"


//Next: References and Further Reading, Previous: Evaluation of B-spline basis function derivatives, Up: Basis Splines
//38.6 Example programs for B-splines

//The following program computes a linear least squares fit to data using cubic B-spline basis functions with uniform breakpoints. The data is generated from the curve y(x) = \cos(x) \exp(-x/10) on [0, 15] with gaussian noise added.


     int smooth_Bspline(QList<float> Xx,  QList<float> Yy, float val)
     {

       int n;
       const size_t ncoeffs = NCOEFFS;
       const size_t nbreak = NBREAK;
       size_t i, j;
       QString temp;
       gsl_bspline_workspace *bw;
       gsl_vector *B;
       double dy;
       gsl_rng *r;
       gsl_vector *c, *w;
       gsl_vector *x, *y;
       gsl_matrix *X, *cov;
       gsl_multifit_linear_workspace *mw;
       double chisq;
       double Rsq;
       double dof;

        n = size_t(Xx.count());
       gsl_rng_env_setup();
       r = gsl_rng_alloc(gsl_rng_default);
        temp.sprintf("%d",Xx.count());
        qDebug(temp.toAscii());
       /* allocate a cubic bspline workspace (k = 4) */
       bw = gsl_bspline_alloc(4, nbreak);
       B = gsl_vector_alloc(ncoeffs);

       x = gsl_vector_alloc(n);
       y = gsl_vector_alloc(n);
       X = gsl_matrix_alloc(n, ncoeffs);
       c = gsl_vector_alloc(ncoeffs);
       w = gsl_vector_alloc(n);
       cov = gsl_matrix_alloc(ncoeffs, ncoeffs);
       mw = gsl_multifit_linear_alloc(n, ncoeffs);

       printf("#m=0,S=0\n");
       /* this is the data to be fitted */
          print2DebugFile("\n\nraw data");
       for (i = 0; i < n; ++i)
         {
           double sigma;
          double xi = Xx.at(i);
          double yi = Yy.at(i);
          /* double sigma;
           double xi = (15.0 / (N - 1)) * i;
           double yi = cos(xi) * exp(-0.1 * xi);

           sigma = 0.1 * yi;
           dy = gsl_ran_gaussian(r, sigma);
           yi += dy;
            */
          gsl_vector_set(x, i, xi);
          gsl_vector_set(y, i, yi);
          temp.sprintf("%d %f %f",i, xi, yi);
          print2DebugFile(temp);
          //temp.sprintf("%f %f\n", xi, yi);
          // qDebug(temp.toAscii());
         }

       /* use uniform breakpoints on [0, 15] */
       gsl_bspline_knots_uniform(Xx.first(),Xx.last(), bw);

       /* construct the fit matrix X */

       for (i = 0; i < n; ++i)
         {
           double xi = gsl_vector_get(x, i);

           /* compute B_j(xi) for all j */
           gsl_bspline_eval(xi, B, bw);

           /* fill in row i of X */
           /*for (j = 0; j < ncoeffs; ++j)
             {
               double Bj = gsl_vector_get(B, j);
               gsl_matrix_set(X, i, j, Bj);
             }
             */
           temp.sprintf("%d %f %f",i, xi, B->data);
             print2DebugFile(temp);
         }

       /* do the fit *//*
       gsl_multifit_linear(X,y, c, cov, &chisq, mw);

       dof = n - ncoeffs;
       Rsq = 1.0 - chisq / gsl_stats_wtss(w->data, 1, y->data, 1, y->size);

       fprintf(stderr, "chisq/dof = %e, Rsq = %f\n", chisq / dof, Rsq);

       /* output the smoothed curve */
       /*
       {
         double xi, yi, yerr;
        temp.sprintf("\n\nsmoothed");
        print2DebugFile(temp);
         printf("#m=1,S=0\n");
         i=0;
         for (xi = Xx.first(); xi < Xx.last(); xi +=  (Xx.last()-Xx.first())/20)
           {
             gsl_bspline_eval(xi, B, bw);
             gsl_multifit_linear_est(B, c, cov, &yi, &yerr);
             temp.sprintf("%d %f %f",i, xi, yi);
             print2DebugFile(temp);
             i++;
             //qDebug(temp.toAscii());
           }
       }
*/
       gsl_rng_free(r);
       gsl_bspline_free(bw);
       gsl_vector_free(B);
       gsl_vector_free(x);
       gsl_vector_free(y);
       gsl_matrix_free(X);
       gsl_vector_free(c);
       gsl_vector_free(w);
       gsl_matrix_free(cov);
       gsl_multifit_linear_free(mw);

       return 0;
     } /* main() */


